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HIGH  ORDER  WENO  SCHEMES  FOR  HAMILTON-JACOBI  EQUATIONS  ON 

TRIANGULAR  MESHES  * 

YONG-TAO  ZHANGt  AND  CHI-WANG  SHU^ 

Abstract.  In  this  paper  we  construct  high  order  weighted  essentially  non-oscillatory  (WENO)  schemes 
for  solving  the  nonlinear  Hamilton- Jacobi  equations  on  two-dimensional  unstructured  meshes.  The  main 
ideas  are  nodal  based  approximations,  the  usage  of  monotone  Hamiltonians  as  building  blocks  on  unstruc¬ 
tured  meshes,  nonlinear  weights  using  smooth  indicators  of  second  and  higher  derivatives,  and  a  strategy  to 
choose  diversified  smaller  stencils  to  make  up  the  bigger  stencil  in  the  WENO  procedure.  Both  third-order 
and  fourth-order  WENO  schemes  using  combinations  of  second-order  approximations  with  nonlinear  weights 
are  constructed.  Extensive  numerical  experiments  are  performed  to  demonstrate  the  stability  and  accuracy 
of  the  methods.  High-order  accuracy  in  smooth  regions,  good  resolution  of  derivative  singularities,  and 
convergence  to  viscosity  solutions  are  observed. 

Key  words,  weighted  essentially  non-oscillatory  schemes,  Hamilton- Jacobi  equations,  high-order  accu¬ 
racy,  unstructured  mesh 

Subject  classification.  Applied  and  Numerical  Mathematics 

1.  Introduction.  In  this  paper,  we  consider  the  numerical  solution  of  Hamilton- Jacobi  (H-J)  equations 

+  =0,  </>(», 0)  =  </io(a;).  (1.1) 

Such  equations  appear  often  in  applications,  such  as  in  optimal  control,  differential  games,  image  processing 
and  computer  vision,  and  geometric  optics.  With  the  popularity  of  level  set  methods  [12]  the  necessity  to 
have  good  algorithms  to  solve  H-J  equations  becomes  even  more  obvious. 

There  have  been  many  papers  in  the  literature  developing  numerical  methods  to  solve  H-J  equations. 
In  certain  sense  H-J  equations  are  easier  to  solve  than  their  conservation  laws  counterparts,  because  the 
solutions  here  are  smoother:  typical  solutions  are  continuous  with  possibly  discontinuous  derivatives. 

For  structured  meshes,  finite  difference  methods  similar  to  those  developed  for  conservation  laws  could 
be  easily  designed.  Thus  we  have  the  earlier  second  order  ENO  schemes  of  Osher  and  Sethian  [12],  higher 
order  essentially  non-oscillatory  (ENO)  schemes  of  Osher  and  Shu  [13],  higher  order  weighted  ENO  (WENO) 
schemes  of  Jiang  and  Peng  [7],  and  central  high  resolution  schemes  of  Lin  and  Tadmor  [10],  among  many 
others. 

However,  for  unstructured  meshes  there  are  concept ional  difficulties  in  designing  schemes  such  as  finite 
volume  schemes  and  finite  element  schemes,  which  are  very  useful  for  conservation  laws.  The  problem  arises 
because  the  H-J  equations  cannot  be  written  in  a  “conservation  form”,  suitable  for  integration  by  parts 
which  is  the  backbone  of  finite  volume  and  finite  element  methods.  Thus  algorithms,  especially  high  order 
algorithms  for  H-J  equations  on  unstructured  meshes  are  relatively  few.  We  have  the  very  good  first  and 
second  order  finite  volume  type  schemes  of  Abgrall  [1]  and  ENO  or  limiter  type  high  order  version  of  Augoula 
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Hampton,  VA  23681-2199,  and  AFOSR  grant  F49620-99- 1-0077. 
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and  Abgrall  [2].  The  monotone  Hamiltonians  developed  in  [1]  are  used  in  this  paper  as  building  blocks.  We 
also  have  the  continuous  finite  element  methods  of  Barth  and  Sethian  [4],  and  the  discontinuous  Galerkin 
methods  of  Hu  and  Shu  [5].  However,  the  formulation  of  the  finite  element  methods  in  [5]  actually  uses  the 
differentiated  version  of  H-J,  which  is  a  system  of  conservation  laws.  It  would  be  desirable  in  many  situations 
to  use  directly  H-J  on  an  unstructured  mesh  and  still  obtain  high  order,  non-oscillatory  numerical  results. 
The  algorithms  developed  in  this  paper  fulfill  this  purpose. 

The  WENO  methodology  adopted  in  this  paper  can  be  traced  back  to  the  earlier  work  for  conservation 
laws  started  in  [11]  for  a  third  order  finite  volume  version  in  one  space  dimension  and  in  [8]  for  third  and 
fifth  order  finite  difference  WENO  schemes  in  multi  space  dimensions  with  a  general  framework  for  the 
design  of  the  smoothness  indicators  and  nonlinear  weights.  The  techniques  most  relevant  to  the  current 
paper  are  the  third  and  fourth  order  finite  volume  WENO  schemes  for  2D  conservation  laws  in  general 
triangulations  [6]  and  the  finite  difference  WENO  schemes  for  Hamilton- Jacobi  equations  [7].  However, 
significant  new  components  of  the  algorithms  have  been  developed  in  this  paper  to  deal  with  the  additional 
difficulty  caused  by  the  “non-conservation”  form  of  the  H-J  equations  and  unstructured  meshes.  These 
include  nodal  based  approximations,  the  usage  of  monotone  Hamiltonians  as  building  blocks  on  unstructured 
meshes,  nonlinear  weights  using  smooth  indicators  of  second  and  higher  derivatives,  and  a  strategy  to  choose 
diversified  smaller  stencils  to  make  up  the  bigger  stencil  in  the  WENO  procedure.  Both  third-order  and 
fourth-order  WENO  schemes  using  combinations  of  second-order  approximations  with  nonlinear  weights 
are  constructed.  Extensive  numerical  experiments  are  performed  to  demonstrate  the  stability  and  accuracy 
of  the  methods.  High-order  accuracy  in  smooth  regions,  good  resolution  of  derivative  singularities,  and 
convergence  to  viscosity  solutions  are  observed. 

The  algorithm  is  developed  in  section  2.  Section  3  contains  numerical  examples  verifying  the  stability, 
convergence  and  accuracy  of  the  algorithms.  Concluding  remarks  are  given  in  section  4. 

2.  The  WENO  algorithm  for  2D  unstructured  meshes.  In  this  section  we  develop  third  and 
fourth  order  WENO  schemes  on  unstructured  meshes  in  2D  for  the  H-J  equations. 

2.1.  The  framework.  We  take  d  =  2  in  (1.1),  and  use  x^y  instead  of  Xi,X2: 

<pt+H{(p^,(py)=0,  <p{x,y,0)  =  (2-1) 


The  equation  (2.1)  is  solved  in  the  domain  which  has  a  triangulation  Th  consisting  of  triangles.  The  nodes 
will  be  named  by  their  indices  0  <  i  <  iV,  with  a  total  of  iV  -h  1  nodes.  For  every  node  i,  we  define  the  + 1 
angular  sectors  To,  •  •  •  ,  Tj^.  meeting  at  the  point  i;  they  are  the  inner  angles  at  node  i  of  the  triangles  having 
i  as  a  vertex.  The  indexing  of  the  angular  sectors  is  ordered  counterclockwise.  is  the  unit  vector  of  the 
half-line  =  T/  flT/+i,  and  9i  is  the  inner  angle  of  sector  T/,  0  <  /  <  kf,  see  Figure  2.1. 

We  will  denote  by  (pi  the  numerical  approximation  to  the  viscosity  solution  of  (2.1)  at  node  i.  (V0)o,  •  •  •  ,  (y4>)ki 
will  respectively  represent  the  numerical  approximation  of  at  node  i  in  each  angular  sector  Tq,  •  •  •  ,  Tj^^. 

An  important  building  block  for  the  algorithms  in  this  paper  is  the  Lax-Friedrichs  type  monotone 
Hamiltonian  for  arbitrary  triangulations  developed  by  Abgrall  in  [1],  which  is  a  generalization  of  the  Lax- 
Friedrichs  monotone  Hamiltonian  for  Cartesian  meshes  in  Osher  and  Shu  [13].  This  monotone  Hamiltonian 
is  given  by 
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Fig.  2.1.  node  %  and  its  angular  sectors. 


where 


A+i  =  tan  M  +  tan  ,  a  =  max{  max  \Hi{u,v)\,  max  \H2{u,v)\}. 


Here  Hi  and  H2  are  the  partial  derivatives  of  H  with  respect  to  4)x  and  (py^  respectively,  or  the  Lipschitz 
constants  of  H  with  respect  to  (px  and  (py,  if  H  is  not  differentiable.  [A^B]  is  the  value  range  for  {(px)h 
and  [C^D]  is  the  value  range  for  (0^)/,  over  0  <  /  <  for  the  local  Lax-Friedrichs  Hamiltonian,  and  over 
0  <  I  <  ki  and  0  <i  <  N  for  global  Lax-Friedrichs  Hamiltonian. 

The  H  in  (2.2)  defines  a  monotone  Hamiltonian.  It  is  Lipschitz  continuous  in  all  arguments  and  is 
consistent  with  i?,  i.e.,  ^(V0,  •••  ,V0)  =  H{V(p).  Hence  if  we  have  high-order  approximations  to  Vcp  at 
node  i  in  every  angular  sector,  the  numerical  Hamiltonian  H  will  be  a  high-order  approximation  to  H. 

The  semi-discrete  scheme  is  thus  given  by: 


-Mt)  +  i?((V</-)o,  •  •  •  ,  (V4>)k,)  =  0  (2.3) 

and  it  is  discretized  in  time  by  the  high  order  nonlinearly  stable  Runge-Kutta  time  discretization  in  [16]. 
If  the  spatial  dimension  reduces  to  one,  the  numerical  Hamiltonian  (2.2)  will  reduce  to  the  local  or  global 
Lax-Friedrichs  Hamiltonian  [13]: 


H{u-,u^)=H 


with  a  =  max  \H\u)\.  If  the  maximum  for  computing  a  is  taken  over  the  range  covered  by  u  and 

A<u<B 

we  get  the  local  Lax-Friedrichs  Hamiltonian;  if  it  is  also  taken  over  all  grid  points,  we  get  the  global 
Lax-Friedrichs  Hamiltonian. 

2.2.  Linear  schemes.  Now  we  discuss  how  to  construct  a  high-order  approximation  for  in  every 
angular  sector  of  every  node.  Let  denote  the  set  of  two-dimensional  polynomials  of  degree  less  than  or 
equal  to  k.  We  use  Lagrange  interpolations  as  follows:  given  a  smooth  function  and  a  triangulation  with 
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triangles  { Aq,  Ai,  . . .  ,  Am}  and  nodes  {0, 1,2,...  ,  N},  we  would  like  to  construct,  for  each  triangle  A^,  a 
polynomial  p{x,y)  G  such  that  p{xi^yi)  =  (t>{xi^yi)^  where  (x/,^//)  are  the  coordinates  of  the  three  nodes 
of  the  triangle  A^  and  a  few  neighboring  nodes.  p{x^y)  would  thus  be  a  (fc  +  l)th-order  approximation  to  0 
on  the  cell  A^. 

Because  kth  degree  polynomial  p{x,y)  has  K  =  degrees  of  freedom,  we  need  to  use  the 

information  of  at  least  K  nodes.  In  addition  to  the  three  nodes  of  the  triangle  A^,  we  may  take  the  other  K—3 
nodes  of  the  neighboring  cells  around  triangle  A^.  We  rename  these  K  nodes  as  Si  =  {Mi,  M2, . . .  ,  M^}, 
Si  is  called  a  big  stencil  for  the  triangle  A^.  Let  (xi^yi)  be  the  barycenter  of  A^.  Define  ^  =  {x  —  Xi)/ hi, 
P—{y  —  yi)/hi^  where  hi  =  y^jA^j  with  jA^j  denoting  the  area  of  the  triangle  A^,  then  we  can  write  p{x,y) 
as: 

k 

Pix,y)  =  '^  XI 

j—0 

Using  the  K  interpolation  conditions: 

p{Mi)  =  mi).  I  =  h2,---,K, 

we  get  d,  K  X  K  linear  system  for  the  K  unknowns  asr-  The  normalized  variables  ry  are  used  to  make  the 
condition  number  of  the  linear  system  independent  of  the  mesh  size. 

It  is  well  known  that  in  two  and  higher  dimensions  this  interpolation  problem  is  not  always  well  defined. 
The  linear  system  can  be  very  ill-conditioned  or  even  singular,  in  such  cases  we  would  have  to  add  more 
nodes  to  the  big  stencil  Si  from  the  neighboring  cells  around  triangle  A^  to  obtain  an  over-determined  linear 
system,  and  then  use  the  least-square  method  to  solve  it. 

After  we  have  obtained  the  approximation  polynomial  p{x,  y)  on  the  triangle  A^,  Vp  will  be  a  fcth-order 
approximation  for  on  A^.  Hence  we  get  the  high-order  approximation  Vp(x/,2//)  to  V0(x/,2//),  for  any 
one  of  the  three  vertices  {xi,yi)  of  the  triangle  A^,  in  the  relevant  angular  sectors. 

2.2.1.  Third-order  linear  schemes.  A  scheme  is  called  linear  if  it  is  linear  when  applied  to  a  linear 
equation  with  constant  coefficients.  We  need  a  third-order  approximation  for  to  construct  a  third-order 
linear  scheme,  hence  we  need  a  cubic  polynomial  interpolation.  A  cubic  polynomial  p^  has  10  degrees  of 
freedom.  We  will  use  some  or  all  of  the  nodes  shown  in  Figure  2.2  to  form  our  big  stencil.  For  our  target 
triangle  Aq,  which  has  three  vertices  k  and  the  barycenter  G,  we  need  to  construct  a  cubic  polynomial 
p^,  then  Vp^  will  be  a  third-order  approximation  for  on  Aq,  and  the  values  of  Vp^  at  points  i,  j  and 
k  will  be  third-order  approximations  for  at  the  angular  sector  Aq  of  nodes  i,  j  and  k.  We  name  the 
nodes  of  the  neighboring  triangles  of  triangle  Aq  as  follows:  nodes  1,  4,  7  are  the  nodes  (other  than  i,  j,  k)  of 
neighbors  of  Aq,  nodes  9,  2,  3,  5,  6,  8  (other  than  1,  4,  7,  i,  j,  k)  are  the  nodes  of  the  neighbors  of  the  three 
neighboring  triangles  of  Aq.  Notice  that  the  points  9,  2,  3,  5,  6,  8  do  not  have  to  be  six  distinct  points.  For 
example  the  points  5  and  6  could  be  the  same  point. 

Our  interpolation  points  are  nodes  i,  j,  fc,  1,  4,  7  and  some  of  the  nodes  taken  from  the  set  E  = 
{9, 2, 3, 5, 6, 8}.  The  algorithm  to  determine  the  big  stencil  So  for  triangle  Aq  is  as  follows. 

Procedure  2.1:  The  choice  of  the  big  stencil  for  the  third-order  scheme. 

1.  Nodes  i,  j,  k,  1,  4,  7  are  always  included  in  So- 

2.  Compute  the  distance  of  point  G  and  every  node  in  the  set  E  =  {9, 2, 3, 5, 6, 8}  respectively.  Sort 
the  nodes  in  E  with  this  distance  from  smallest  to  biggest,  and  denote  the  sorted  node  set  by  E. 
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Fig.  2.2.  The  nodes  used  for  the  big  stencil  of  the  third-order  scheme. 


3.  From  the  sorted  node  set  E,  add  the  first  4  nodes  to  Sq.  Use  the  10  nodes  in  So  as  interpolation 
points  to  get  the  10  x  10  interpolation  coefficient  matrix  A, 
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where  =  (x/  -  XG)/hG,  rji  =  {yi  -  yG)/hG,  hG  =  \/|Ao|,  and  (xi.yi)  are  the  coordinates  of  the 
nodes  in  5o,  (xG^yo)  is  the  coordinate  of  the  point  G. 

4.  Compute  the  reciprocal  condition  number  c  of  A.  This  is  provided  by  most  linear  solvers.  If  c>  6 
for  some  threshold  5,  we  have  obtained  the  final  stencil  Sq.  Otherwise,  add  the  fifth  node  in  E  to 
So-  Use  the  11  nodes  in  So  as  interpolation  points  to  get  the  11  x  10  least  square  interpolation 
coefficient  matrix  A.  Judge  the  reciprocal  condition  number  c  again.  Continue  in  doing  this  until 
c>  d  is  satisfied.  In  our  computation,  we  have  taken  d  =  10“^  as  a  good  threshold  after  extensive 
numerical  experiments.  Notice  that,  since  we  have  normalized  the  coordinates,  this  threshold  does 
not  change  when  the  mesh  is  scaled  uniformly  in  all  directions.  For  all  the  triangulations  we  have 
tested,  at  most  12  nodes  are  needed  in  So  to  reach  the  condition  c>  6.  ■ 

Remark  2.1.  The  reason  that  we  always  include  the  nodes  i,  j,  fc,  1,  4,  7  in  So  is  to  make  the  target  cell 
Ao  sufficiently  “central”  in  the  stencil,  thus  avoiding  So  to  be  seriously  downwind  biased  which  could  lead 
to  linear  instability.  Linear  instability  is  indeed  observed  in  our  numerical  experiments  when  1,  4,  7  are  not 
forced  to  be  in  5o.  ■ 

We  now  have  obtained  the  big  stencil  So  and  its  associated  cubic  polynomial  p^.  For  each  node  (xi^yi) 
in  Ao,  'Vp^{xi,yi)  is  a  third-order  approximation  to  V(p{xi,yi).  In  order  to  construct  a  high-order  WENO 
scheme,  an  important  step  is  to  obtain  a  high-order  approximation  using  a  linear  combination  of  lower  order 
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approximations.  We  will  use  a  linear  combination  of  second-order  approximations  to  get  the  same  third-order 
approximation  to  V0(x/,2//)  as  he.,  we  require 

^P^{xuyi)  =  Y,ls,x^Ps{xuyi),  -^P^ixi,yi)  =  Y,^7s,y-^Ps(xi,yi)  (2.5) 

where  Ps  are  quadratic  interpolation  polynomials,  and  ^s,x  and  ^s,y  are  the  linear  weights  for  the  x-directional 
derivative  and  the  ^-directional  derivative  respectively,  for  s  =  1,  •  •  •  ,q.  The  linear  weights  are  constants 
depending  only  on  the  local  geometry  of  the  mesh.  The  equalities  in  (2.5)  should  hold  for  any  choices  of  the 
function  0. 

Notice  that  to  get  a  second-order  approximation  for  the  derivatives  V0(x/,2//),  we  need  a  quadratic 
interpolation  polynomial.  According  to  the  argument  in  [6],  the  cubic  polynomial  p^{x,y)  has  four  more 
degrees  of  freedom  than  each  quadratic  polynomial  Ps{x,y)^  namely  ^x^y^xy^^ ^y^ .  For  the  six  degrees  of 
freedom  l^x^y^x‘^ ^xy^y‘^ ^  if  we  take  cj)  =  l^cj)  =  x^cj)  =  y^cj)  =  x^^^cj)  =  xy  and  0  =  2/^,  the  equalities  in  (2.5) 
will  hold  for  all  these  cases  under  only  one  constraint  each  on  ^s,x  and  7^,^,  namely  Y^l^i"fs,x  =  1  and 
Yll=i  ls,y  =  I5  because  p^  and  Ps  all  reproduce  these  functions  exactly.  Hence  we  should  only  need  ^  >  5. 
We  take  ^  =  5  in  our  scheme. 

We  now  need  ^  =  5  small  stencils  F^,  s  =  1,  •  •  •  ,  5  for  the  target  triangle  Aq,  satisfying  So  =  Us=i 
and  every  quadratic  polynomial  Ps  is  associated  with  a  small  stencil  F^.  In  our  third-order  scheme,  the  small 
stencils  will  be  the  same  for  both  directions  x,  y  and  all  three  nodes  i^j^k  in  Aq.  However  the  linear  weights 
ls,x^ls,y  can  be  different  for  different  nodes  i^j^k  and  different  directions  x^y.  Because  each  quadratic 
polynomial  has  six  degrees  of  freedom,  the  number  of  nodes  in  F^  must  be  at  least  six.  To  build  a  small 
stencil  F^,  we  start  from  several  candidates  T^J'\r  =  1, 2,  •  •  •  ,77,5.  These  candidates  are  constructed  by  first 
taking  a  point  as  the  “center” ,  then  finding  at  least  six  nodes  from  So  which  have  the  shortest  distances 
from  and  can  generate  the  interpolation  coefficient  matrix  with  a  good  condition  number,  using  the 
method  of  Procedure  2.1.  We  then  choose  the  best  F^  among  Ti^\r  =  !,•••  , for  every  s  =  1,  •  •  •  ,5.  Here 
“best”  means  that  by  using  this  group  of  small  stencils,  the  linear  weights  'ys,x:ls,y:  s  =  1,  •  •  •  ,  5  for  all  three 
nodes  i,  j,  k  are  either  all  positive  or  have  the  smallest  possible  negative  values  in  magnitude.  The  details  of 
the  algorithm  is  described  in  the  following  procedure. 

Procedure  2.2:  The  third-order  linear  scheme. 

For  every  triangle  A/,  /  =  !,•••  ,  iV,  do  steps  1  ~  5: 

1.  Follow  Procedure  2.1  to  obtain  the  big  stencil  5/  for  A/. 

2.  For  s  =  1,  •  •  •  ,  5,  find  the  set  Wg  =  r  =  1, 2,  •  •  •  ,  n^},  which  are  the  candidate  small  stencils 
for  the  s-th  small  stencil.  We  use  the  following  method  to  find  the  F^^^  in  Wgi  first,  nodes  i^j^k 
are  always  included  in  every  F^^^ ;  then  we  take  a  point  A^g^  as  the  center  of  F^^^ ,  detailed  below, 
and  using  the  idea  of  Procedure  2.1,  we  find  at  least  3  additional  nodes  other  than  i,  j,  from  Si 
which  satisfy  the  following  two  conditions:  1.)  they  have  the  shortest  distances  from  A^f^;  and 
2.)  taking  them  and  the  nodes  i^j^k  as  the  interpolation  points,  we  will  obtain  the  interpolation 
coefficient  matrix  A  with  a  good  condition  number,  namely  the  reciprocal  condition  number  c  of 
A  satisfies  c  >  6  with  the  same  threshold  6  =  10“^.  For  the  triangulations  we  have  tested,  at 
most  8  nodes  are  used  to  reach  this  threshold  value.  Finally,  the  center  of  the  candidate  stencils 
A^J'\r  =  1,  •  •  •  , s  =  1,  •  •  •  ,  5  are  taken  from  the  nodes  around  A/  (see  Figure  2.2)  as  follows: 

•  A^^^  =  point  G,  ni  =  1; 

•  A^^  =  node  1,  A^^  =  node  2,  A^^  =  node  9,  77,2  =  3; 
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•  =  node  4,  =  node  3,  ^3^^  =  node  5,  77,3  =  3; 

•  =  node  7,  A^"^  =  node  6,  A^"^  =  node  8,  77,4  =  3; 

•  =  nodes  2,9,3, 5,6,8  and  the  middle  points  of  nodes  2  and  3,  5  and  6,  9  and  8. 
77,5  =  9. 

3.  By  taking  one  small  stencil  from  each  Ws,  s  =  1,  •  •  •  ,  5  to  form  a  group,  we  obtain  77,1  x  77,2  x 
•  ••  X  77,5  groups  of  small  stencils.  We  eliminate  the  groups  which  contain  the  same  small  stencils, 
and  also  eliminate  the  groups  which  do  not  satisfy  the  condition 

U  =  Si 

s—1 

According  to  every  group  =  !,•••  ,5}  of  small  stencils,  we  have  5  quadratic  polynomials 

We  evaluate  and  at  points  to  obtain  second-order  approximation 

values  for  at  the  three  vertices  of  the  triangle  A/.  We  remark  that  for  practical  implementation, 
we  do  not  use  the  polynomial  itself,  but  compute  a  series  of  constants  which  depend  on  the 

local  geometry  only,  such  that: 


dx 


/=1 


(2.6) 


where  every  constant  ai  corresponds  to  one  node  in  the  stencil  F^^^^  and  m  is  the  total  number  of 
nodes  in  For  every  vertex  (xn^Vn)  of  triangle  A/,  we  obtain  a  series  of  such  constants.  And 

for  the  y  directional  partial  derivative,  we  compute  the  corresponding  constants  too. 

4.  For  every  group  s  =  1,  •  •  •  ,  5},  we  form  linear  systems  and  solve  them  to  get  a  series  of  linear 

weights  and  satisfying  the  equalities  (2.5),  for  the  three  vertices  i,  j,  k.  Using  the  previous 
argument  for  combining  low-order  approximations  to  get  a  high-order  approximation,  we  form  the 
linear  system  for  at  a  vertex  {in^Pn)  as  follows  (note  that  we  use  normalized  variables):  take 
4>  —  respectively,  the  equalities  are: 


(2.7) 


where  p^f^^  is  the  quadratic  interpolation  polynomial  for  0,  using  stencil  F^^^^  Again,  in  practical 

(r  ) 

implementation,  we  will  not  use  ph  itself,  instead  we  use  the  constants  computed  in  the  last 
step  and  equation  (2.6)  to  compute  the  approximation  for  the  derivatives  of  0.  Together  with  the 
requirement 


5 


5=1 


(2.8) 


we  obtain  a  5  x  5  linear  system  for  For  the  same  argument  can  be  applied.  Note  that 

we  need  to  compute  the  reciprocal  condition  number  c  for  every  linear  system  again.  If  c  >  5  for  the 
same  threshold  6  =  10“^,  we  will  accept  this  group  of  stencils  as  one  of  the  remaining  candidates. 
Otherwise,  the  linear  system  is  considered  to  be  ill-conditioned  and  its  corresponding  group  of  small 
stencils  s  =  !,•••  ,  5}  is  eliminated  from  further  consideration. 

5.  For  each  of  the  remaining  group  A/  =  s  =  1,  •  •  •  ,  5},  find  the  minimum  value  7/  of  all  these 

linear  weights  of  the  three  vertices  i,  j,  k.  Then  find  the  group  of  small  stencils  whose  7/ 
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Fig.  2.3.  The  nodes  used  for  the  big  stencil  of  the  fourth-order  scheme. 


is  the  biggest  among  all  {7/}^i,  and  take  this  group  as  our  final  5  small  stencils  for  triangle  A/. 
Denote  them  by  s  =  1,  •  •  •  ,  5.  For  every  final  small  stencil  F^,  s  =  1, 2,  •  •  •  ,  5,  we  store  the  index 
numbers  of  the  nodes  in  F^,  the  constants  in  the  linear  combinations  of  node  values  to  approximate 
values  of  at  points  i,  fc,  and  the  linear  weights  ^s,x:  ls,y  oi  the  three  points  i,  k. 

6.  Now  we  have  set  up  the  necessary  constants  which  only  depend  on  the  mesh  for  all  triangles.  To  form 
the  final  linear  scheme,  we  compute  the  third-order  approximations  (V0)o,  •  •  •  ?  (V0)a;^  for  all  mesh 
nodes  1,  by  the  linear  combinations  of  second-order  approximations,  using  the  prestored  constants 
and  linear  weights.  Then  we  can  form  the  scheme  (2.3).  ■ 

Remark  2.2.  The  strategy  described  above  of  finding  second  order  smaller  stencils  to  make  up  the  third 
order  big  stencil  has  been  reached  after  our  extensive  numerical  experiments.  They  are  found  to  be  robust  to 
different  triangulations  and  equations.  Of  course  we  are  not  claiming  that  this  procedure  will  not  fail.  We  can 
only  claim  that  it  behaves  nicely  in  all  our  numerical  tests.  The  computational  cost  of  this  procedure  is  quite 
high,  as  many  choices  and  comparisons  are  needed.  Fortunately,  at  least  for  problems  without  adaptive  mesh 
refinements,  this  procedure  is  done  only  once  at  the  beginning,  and  does  not  need  to  be  repeated  during  time 
marching.  Because  of  the  prestored  constant  coefficients  for  computing  the  the  approximation  to  derivatives, 
the  time  evolution  part  of  the  scheme  is  very  efficient.  ■ 

2.2.2.  Fourth-order  linear  schemes.  To  construct  a  fourth-order  linear  scheme,  we  need  a  fourth- 
order  approximation  to  V0.  Hence  a  fourth  degree  interpolation  polynomial  p^{x,y)  is  required,  which  has 
15  degrees  of  freedom.  We  will  use  part  or  all  of  the  nodes  shown  in  Figure  2.3  to  construct  the  big  stencil 
of  the  fourth-order  scheme. 

Comparing  Figure  2.3  with  Figure  2.2,  we  can  see  that,  in  order  to  get  the  set  of  candidate  nodes  for 
the  big  stencil  of  the  fourth-order  scheme,  we  have  added  nodes  10, 11, 12, 13, 14, 15  (Figure  2.3)  to  the  set 
for  the  third  order  scheme  (Figure  2.2).  Notice  again  that  some  of  the  added  nodes  may  not  be  distinct. 

The  procedure  to  determine  the  big  stencil  of  the  fourth-order  scheme  is  slightly  different  from  that  for 
the  third-order  scheme.  This  difference  is  mainly  due  to  implementation  efficiency  considerations  and  the 
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more  stringent  stability  requirement  for  the  higher  order  scheme.  We  will  not  use  the  distances  between  the 
barycenter  point  and  the  nodes  to  sort  our  candidate  nodes  as  in  Procedure  2.1.  Instead,  we  give  an  ordering 
to  the  candidate  nodes.  This  ordering  is  given  so  that,  when  the  nodes  are  chosen  sequentially  from  it  to 
form  the  big  stencil  5o,  the  target  triangle  Aq  remains  central  to  avoid  serious  downwind  bias  which  could 
lead  to  linear  instability.  Referring  to  Figure  2.3,  our  interpolation  points  for  the  polynomial  include 
nodes  i,  j,  k  and  the  nodes  taken  from  the  sorted  set:  W  =  {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15}.  The 
detailed  procedure  to  determine  the  big  stencil  So  for  the  target  triangle  Aq  is  given  below. 


Procedure  2.3:  The  big  stencil  of  the  fourth-order  scheme. 

1.  Find  15  different  nodes  (other  than  the  nodes  i,j,  fc)  around  triangle  Aq,  and  give  the  order  for 
them  as  in  Figure  2.3.  Naming  them  using  their  ordering  number,  we  get  the  sorted  node  set: 
W  =  {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15}. 

2.  To  start  with,  we  take  So  =  {i,  j,  fc,  1, 2,3,4, 5,6, 7, 8,9, 10, 11, 12}.  Use  this  stencil  So  to  form  the 
15  X  15  interpolation  coefficient  matrix  A  as  in  Procedure  2.1. 

3.  Compute  the  reciprocal  condition  number  c  of  If  c  >  5  for  the  same  threshold  S  =  10“^,  we  have 
obtained  the  final  big  stencil  So-  Otherwise,  we  add  the  next  node  (i.e.  node  13)  in  W  to  So-  Then 
we  obtain  the  16  x  15  interpolation  coefficient  matrix  A  associated  with  this  5o,  an  over-determined 
system  which  can  be  solved  by  the  least  square  procedure.  Compute  the  reciprocal  condition  number 
of  A  and  check  whether  c>  6  again.  Repeat  this  if  necessary  until  c  >  5  is  satisfied.  In  our  numerical 
tests,  at  most  16  nodes  are  needed  in  So  to  reach  our  threshold  value  in  all  the  cases.  ■ 

We  emphasize  again  that  the  ordering  for  the  nodes  as  in  Figure  2.3  is  important  to  guarantee  the  chosen 
stencil  to  yield  to  linear  stability  of  the  scheme. 

Using  the  big  stencil,  we  obtain  the  fourth-order  approximation  Vp^  for  V0.  The  key  step  in  building 
our  fourth-order  WENO  scheme  is  to  get  a  fourth-order  approximation  for  based  on  lower  order  approx¬ 
imations.  We  still  would  like  to  construct  several  second-order  approximations  whose  weighted  average  will 
give  the  same  result  as  Vp^  at  each  angular  sector  of  every  node,  i.e.. 


where  Ps,x:Ps,y  sire  the  quadratic  interpolation  polynomials,  and  sire  the  linear  weights. 

Remark  2.3.  The  reason  that  we  still  use  second-order  approximations,  not  third-order  ones,  as  the  lower 
order  approximations  is  that  each  second-order  approximation  needs  fewer  nodes  in  the  big  stencil,  thus  it  is 
easier  to  make  the  small  stencils  to  be  sufficiently  diversified,  which  is  important  for  the  WENO  technique  to 
function  well  near  discontinuous  derivatives  in  the  solutions.  We  have  encountered  difficulties  in  obtaining 
non-oscillatory  results  converging  to  the  correct  viscosity  solutions  for  some  of  the  more  demanding  test 
cases  when  third  order  building  blocks  are  used  instead  of  the  second  order  ones.  ■ 

We  now  determine  the  value  of  which  is  the  number  of  small  stencils,  in  the  equalities  (2.9).  Because 
has  nine  more  degrees  of  freedom  than  the  quadratic  polynomials  Ps,x^Ps,y^  be.,  xy^^  ^x^  ^x^y^  x^y^^ 

xy^^y^^  according  to  our  previous  argument,  we  would  need  q  >  10.  We  take  ^  =  10  in  our  fourth-order 
scheme.  The  procedure  to  find  all  the  small  stencils  is  described  below. 


Procedure  2.4:  The  fourth-order  linear  scheme. 

For  every  triangle  A/,  /  =  1,  •  •  •  ,  iV,  do  steps  1  ~  3: 

1.  Follow  Procedure  2.3  to  obtain  the  big  stencil  Si  for  A/. 
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2.  For  each  of  the  six  approximations  (py  at  vertices  i,  k  of  A/),  find  10  small  stencils  respectively. 
Let  {Fg^a;?  -s  =  1,  •  •  •  ,  10}  denote  the  small  stencils  to  approximate  (px  sit  vertex  i.  To  find  them,  we 
first  determine  10  preliminary  small  stencils  {F^,  s  =  1,  •  •  •  ,  10}  as  follows  (see  Figure  2.3): 

(1)  Include  nodes  {i,  j,  fc,  1,4, 7}  in  F^. 

(2)  Include  nodes  {i,  j,  k,  1, 2, 9}  in  F^. 

(3)  Include  nodes  {i,  j,  fc,4,3, 5}  in  Fg. 

(4)  Include  nodes  {i,  j,  fc,  7, 8, 6}  in  F^. 

(5)  Include  nodes  {i,  j,  k}  in  all  F^,  s  =  5, 6,  •  •  •  ,  10. 

(6)  Take  points  As,  s  =  1,  •  •  •  ,  10  as  the  center  of  {r^}^2zi5  where  Ai  =  point  G,  A2  =  node  1,  As  = 
node  4,  A4  =  node  7,  A^  =  node  10,  Aq  =  node  13,  ^7  =  node  11,  As  =  node  14,  Ag  =  node 
12,  and  Aiq  =  node  15. 

(7)  Obtain  6  nodes  for  each  of  {r^}^?.!-  There  are  already  6  points  in  For  {r^}g2.5,  we 

add  3  nodes  (other  than  i,j,  k)  from  Si  to  each  of  them,  which  have  the  shortest  distances  from 
As.  Then  use  F^  to  get  the  interpolation  coefficient  matrix,  and  judge  the  reciprocal  condition 
number  c  of  the  matrix.  If  c  reaches  our  threshold  value  6  =  10“^,  then  F^  is  found.  Otherwise 
add  one  different  node  from  5/  to  F^,  which  has  the  shortest  distance  from  As.  Continue  this 
until  c>  6  is  satisfied. 

(8)  All  of  the  6  approximations  (to  4)x,  (py  at  vertices  i,j,  k  of  A/)  have  the  common  10  preliminary 
small  stencils  {F^,  s  =  1,  •  •  •  ,  10}. 

Now  we  come  to  determine  the  {Ts^x:  =  1?  *  *  *  5 10}-  go^il  is  that  the  coefficient  matrix  A  of 

the  10  X  10  linear  system,  which  is  obtained  from  and  used  to  compute  the  linear  weights, 

has  a  good  condition  number.  For  s  =  l,2,---,10,  we  perform  the  following  steps: 

(1)  Taking  the  original  small  stencil  F^  as  the  interpolation  stencil,  compute  the  constants  {a/}^i, 
which  depend  on  the  mesh  only  and  are  the  coefficients  in  the  linear  combination  of  function 
values  at  the  nodes  in  F^  to  get  the  second-order  approximation  to  (px  at  vertex  i. 

(2)  Form  the  s-th  column  of  the  matrix  A.  Let  art  be  the  elements  of  the  matrix  A.  Take  the 

9  functions  {0(^^}}2z2  =  respectively,  and  let  p^J^l  be  the 

second-order  interpolation  polynomial  for  (p^'^\r  =  2,  •  •  •  ,10,  then  ai^  =  1,  and 

are  =  r  =  2,...,10  (2.10) 

where  {^i,r}i)  is  the  coordinate  of  vertex  i  (again  note  that  we  use  normalized  variables).  We  use 
the  constants  computed  at  last  step  to  implement  this,  and  will  not  need  to  use  the  polynomials 
themselves. 

(3)  Now  A  has  s  columns  and  is  a  10  x  s  matrix.  Compute  its  reciprocal  condition  number  c.  If 
c>  d,  take  the  current  small  stencil  as  our  s-th  small  stencil  Ts^x-  Otherwise,  change  one  node 
which  is  in  the  current  small  stencil  and  farthest  from  the  center  point  As,  to  another  node 
from  Si  which  is  not  in  the  current  small  stencil  but  is  nearest  to  As.  Now  we  get  a  new  small 
stencil.  This  part  can  be  repeated  if  c  >  5  is  still  not  satisfied.  Then  the  current  small  stencil 
is  our  s-th  small  stencil. 

For  the  small  stencils  to  approximate  the  y-directional  derivative  at  vertex  i  and  x,y-directional 
derivatives  at  other  2  vertices  j,  k,  we  use  a  similar  procedure. 

3.  Now  we  have  obtained  the  coefficient  matrix  A  with  a  good  condition  number.  Along  with  the  right 
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r  =  2,---  ,10, 


(2.11) 


hand  vector  b,  whose  components  are  hi  =  l  and 

K  = 

we  obtain  the  10  x  10  linear  system 

I7  =  b.  (2.12) 

We  solve  it  to  get  the  linear  weights  7  =  (71, a;?  •  •  •  :7io,x)^’  We  use  similar  method  to  get  {7s,y}52ii- 
Then  for  each  of  small  stencils,  we  store  the  index  numbers  of  the  nodes  in  it,  the  constants  in 
the  linear  combinations  of  node  values  to  approximate  values  of  the  derivatives  of  0,  and  the  linear 
weights. 

4.  Now  we  have  set  up  necessary  constants  which  only  depend  on  the  mesh  for  all  triangles.  To  form  the 
final  linear  scheme,  we  compute  the  fourth-order  approximation  (V(/>)o,  •  •  •  for  all  the  mesh 

nodes  1,  by  the  linear  combinations  of  second-order  approximations,  using  the  prestored  constants 
and  linear  weights.  We  then  form  the  scheme  (2.3).  ■ 

Remark  2.4.  We  notice  that  in  the  fourth-order  scheme,  although  the  big  stencil  is  the  same  to  approximate 
-^(j)  and  at  all  three  vertices  i^j^k  of  the  target  cell  Aq,  the  small  stencils  can  be  different  for  the  x- 
direction  and  y-direction  derivatives  at  one  vertex,  and  also  different  at  the  three  vertices  i,  j,  in  the  same 
cell  Aq.  This  strategy  is  different  from  that  for  the  third-order  scheme  described  before.  The  reason  for  this 
difference  is  that  we  need  more  small  stencils  than  in  the  third-order  case  and  it  is  difficult  to  find  a  group 
of  common  small  stencils  for  all  the  six  cases  (three  vertices,  2  derivatives  each).  ■ 

Remark  2.5.  Again,  the  procedure  described  above  of  finding  second  order  smaller  stencils  to  make  up 
the  fourth  order  big  stencil  has  been  reached  after  extensive  numerical  experiments.  They  are  found  to  be 
robust  to  different  triangulations  and  equations.  Of  course  we  are  not  claiming  that  this  procedure  will 
not  fail.  We  can  only  claim  that  it  behaves  nicely  in  all  our  numerical  tests.  The  computational  cost  of 
this  procedure  is  again  quite  high,  as  many  choices  and  comparisons  are  needed.  Fortunately,  at  least  for 
problems  without  adaptive  mesh  refinements,  this  procedure  is  done  only  once  at  the  beginning,  and  does 
not  need  to  be  repeated  during  time  marching.  Because  of  the  prestored  constant  coefficients  for  computing 
the  the  approximation  to  derivatives,  the  time  evolution  part  of  the  scheme  is  again  very  efficient.  ■ 


2.3.  WENO  schemes.  In  this  section,  we  construct  the  WENO  schemes  based  on  non-linear  weights. 
The  resulting  schemes  will  be  suitable  to  compute  the  H-J  equations  whose  solutions  are  not  smooth. 


2.3.1 .  WENO  approximation.  We  only  discuss  the  case  of  WENO  approximation  for  the  x-directional 
derivative  at  vertex  i  of  the  target  cell  A/.  Other  cases  are  similar.  In  order  to  compute  the  non-linear  weights, 
we  need  to  compute  the  smoothness  indicators  first. 

For  a  polynomial  j9(x,  y)  defined  on  the  target  cell  Aq  with  degree  up  to  fc,  we  take  the  smoothness 
indicator  (3  as: 


/3  = 


{D°‘p{x,  y)f  dxdy 


(2.13) 


where  a  is  a  multi-index  and  D  is  the  derivative  operator.  The  smoothness  indicator  measures  how  smooth 
the  function  p  is  on  the  triangle  Aq:  the  smaller  the  smoothness  indicator,  the  smoother  the  function  p  is  on 
Aq.  The  scaling  factor  in  front  of  the  derivatives  renders  the  smoothness  indicator  self-similar  and  invariant 
under  uniform  scaling  of  the  mesh  in  all  directions. 
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Remark  2.6.  The  definition  of  the  smoothness  indicator  in  equation  (2.13)  is  different  from  that  used  in 
the  WENO  schemes  for  conservation  laws  [8].  In  equation  (2.13),  the  range  of  summation  is  2  <  |a|  <  fc, 
but  in  [8]  for  conservation  laws,  it  is  1  <  |a|  <  fc.  An  intuitive  reason  is  that  H-J  equations  are  in  some  sense 
the  conservation  laws  “integrated  once” ,  hence  smooth  indicators  for  the  former  should  involve  derivatives 
one  order  higher  than  those  for  the  latter.  A  similar  strategy  is  also  used  in  ENO  schemes  for  H-J  equations 
in  [13].  We  have  found  through  numerical  experiments,  that  if  we  still  take  1  <  |a|  <  k  to  compute  the 
smoothness  indicators,  we  cannot  obtain  the  correct  viscosity  solutions  for  some  non-convex  Hamilton- Jacobi 
equations.  ■ 

Now  we  define  the  non-linear  weights  as: 


Ij 


(2.14) 


where  7j  is  the  jth  linear  weight  (e.g.  the  js,x  in  our  linear  schemes),  is  the  smoothness  indicator  for 
the  jth  interpolation  polynomial  Pj{x,y)  (e.g.  the  Ps  in  equation  (2.5)  for  the  third-order  case  and  the  Ps,x 
in  equation  (2.9)  for  the  fourth-order  case)  associated  with  the  jth  small  stencil,  and  £  is  a  small  positive 
number  to  avoid  the  denominator  becoming  0.  We  take  e  =  10“ ®  for  all  the  computations  in  this  paper. 
The  final  WENO  approximation  for  the  x-directional  derivative  at  vertex  i  of  target  cell  A/  is  given  by 


d 

,)i  =  Y^Wj—pj{xi,yi) 


(2.15) 


where  {xi^yi)  are  the  coordinates  of  vertex  i,  ^  =  5  for  the  third-order  schemes  and  ^  =  10  for  the  fourth-order 
schemes. 

In  our  WENO  schemes,  the  linear  weights  depend  on  the  local  geometry  of  the  mesh  and  can 

be  negative.  If  min(7i,  •  •  •  ,7^)  <  0,  we  adopt  the  splitting  technique  of  treating  negative  weights  in  WENO 
schemes  developed  by  Shi,  Hu  and  Shu  [15]:  first  we  split  the  linear  weights  into  two  groups: 

%^  =  ^(7j +3|7j|),  i  =  (2.16) 

then  we  scale  them  by 

i  =  1,  •  •  •  ,  9-  (2.17) 

1=1 

Now  we  compute  the  nonlinear  weights  (2.14)  for  the  positive  and  negative  groups  7^  separately,  denoted 
by  based  on  the  same  smoothness  indicator  pj.  Then  we  compute  the  WENO  approximations  {px)^ 
separately  by  (2.15),  using  and  form  the  final  WENO  approximation  by 


-O-  {(l)x)i  ’ 


(2.18) 


The  key  idea  of  this  decomposition  is  to  make  sure  that  every  stencil  has  a  significant  representation  in 
both  the  positive  and  the  negative  weight  groups.  Within  each  group,  the  WENO  idea  of  redistributing  the 
weights  subject  to  a  fixed  sum  according  to  the  smoothness  of  the  approximation  is  still  followed  as  before. 
See  [15]  for  more  details. 

Again,  we  remark  that  the  smoothness  indicator  (2.13)  is  a  quadratic  function  of  function  values  on 
nodes  of  the  small  stencil,  so  in  a  practical  implementation,  to  compute  the  smoothness  indicator  pj  for  the 
j-th  small  stencil  by  equation  (2.13),  we  do  not  need  to  use  the  interpolation  polynomial  itself,  instead  we 
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use  a  series  of  constants  {a^t,r  =  !,•••  =  !,•••  ,m},  which  can  be  precomputed  and  they  depend  on 

the  mesh  only,  such  that 

m  t 

f3j  =Y,<PtCE^rt<Pr),  (2.19) 

where  m  is  the  total  number  of  nodes  in  the  j-th  small  stencil.  These  constants  for  all  smoothness  indicators 
should  be  precomputed  and  stored  once  the  mesh  is  generated. 

2.3.2.  Algorithm  flowchart.  We  summarize  the  algorithm  for  the  third-order  and  the  fourth-order 
WENO  schemes  as  follows: 

Procedure  2.5:  The  third-  and  fourth-order  WENO  schemes. 

1.  Generate  a  triangular  mesh. 

2.  Compute  and  store  all  constants  which  only  depend  on  the  mesh  and  the  accuracy  order  of  the 
scheme.  These  constants  include  the  node  index  numbers  of  each  small  stencil,  the  coefficients  in 
the  linear  combinations  of  function  values  on  nodes  of  small  stencil  to  approximate  the  derivative 
values  and  the  linear  weights,  following  the  Procedure  2.2  for  the  third-order  case  and  the  Procedure 
2.4  for  the  fourth-order  case,  and  the  constants  for  computing  smoothness  indicators  in  equation 
(2.19). 

3.  Using  the  prestored  constants,  for  each  angular  sector  of  every  node  i,  compute  the  low-order  ap¬ 
proximations  for  V(/>  and  nonlinear  weights,  then  compute  the  third  or  fourth-order  WENO  approx¬ 
imations  (2.15)  or  (2.18).  Form  the  scheme  (2.3).  Use  high-order  TVD  Runge-Kutta  time  stepping 
[16]  to  evolve  in  time. 

3.  Numerical  examples.  In  this  section,  we  apply  the  WENO  schemes  developed  in  the  previous 
section  to  a  set  of  one  and  two  dimensional  problems.  The  CEL  number  is  taken  as  0.5  in  all  the  cases, 
except  for  the  accuracy  test  where  it  is  taken  to  be  smaller  if  necessary  to  guarantee  that  spatial  errors 
dominate.  For  the  temporal  discretization,  we  use  the  third-order  TVD  Runge-Kutta  scheme  of  Shu  and 
Osher  in  [16]. 

We  have  not  described  the  details  of  the  one  dimensional  algorithm  in  the  previous  section.  The  algorithm 
in  ID  is  just  the  same  2D  algorithm  with  only  two  “angular  sectors”.  The  WENO  interpolation  follows  along 
the  lines  of  [8]  and  [3].  In  all  the  one-dimensional  numerical  examples,  we  use  non-uniform  meshes.  These 
non-uniform  meshes  are  obtained  by  randomly  shifting  the  cell  boundaries  in  a  uniform  mesh  in  the  range 
[— O.l/i,  O.l/i]  where  h  is  the  uniform  mesh  size.  When  we  perform  the  accuracy  test,  the  refinement  of  the 
meshes  is  achieved  by  cutting  each  cell  into  two  smaller  equal  sized  ones. 

Example  3.1.  One-dimension  linear  equation: 

f  +  00.  =  0,  0  <  X  <  27r, 

1  0(x,O)  =  sin(x) 

with  periodic  boundary  conditions.  We  remark  that  even  though  this  equation  is  also  a  conservation  law, 
the  method  used  here  is  different  from  the  WENO  schemes  for  conservation  laws  in  [8]. 

We  use  third,  fifth  and  seventh  order  linear  and  WENO  schemes  to  compute  the  problem  to  t  =  2.0. 
The  errors  and  the  numerical  orders  of  accuracy  are  listed  in  Table  3.1.  We  can  see  that  the  correct  orders  of 
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Table  3.1 

Accuracy  for  ID  linear  equation,  Linear  and  WENO  Schemes,  t=2 


Linear 

3rd  order 

5th  order 

7th  order 

N 

LOO  QYYQY 

order 

LOO 

order 

LOO 

order 

10 

4.38E-02 

4.94E-03 

1.23E-03 

— 

20 

5.85E-03 

2.91 

1.68E-04 

4.88 

1.06E-05 

6.86 

40 

7.44E-04 

2.97 

5.56E-06 

4.92 

8.63E-08 

6.94 

80 

9.41E-05 

2.98 

1.79E-07 

4.95 

6.91E-10 

6.96 

160 

1.18E-05 

2.99 

5.67E-09 

4.98 

5.46E-12 

6.98 

WENO 

3rd  order 

5th  order 

7th  order 

N 

LOO 

order 

LOO 

order 

LOO 

order 

10 

1.55E-01 

— 

2.00E-02 

— 

3.07E-03 

— 

20 

4.64E-02 

1.74 

8.64E-04 

4.53 

3.45E-05 

6.48 

40 

1.22E-02 

1.93 

2.76E-05 

4.97 

4.24E-07 

6.35 

80 

1.66E-03 

2.88 

8.75E-07 

4.98 

4.04E-09 

6.71 

160 

5.38E-05 

4.95 

2.33E-08 

5.23 

1.04E-11 

8.61 

accuracy  are  obtained  by  both  the  linear  and  WENO  schemes.  In  all  the  accuracy  tests,  we  only  list  results 
for  errors.  Those  for  or  L?  errors,  which  follow  the  same  patterns,  are  omitted  to  save  space. 


Example  3.2.  One-dimensional  linear  equation: 


(pf  “h  (f)j^  —  0,  1  ^  X  <C  1, 

(p{x,0)  =g(x-0.5) 
with  periodic  boundary  conditions.  Where 

"2cos(^)- v^. 


.  .  .\/3  9  27r.  .  _  j 

9ix)  =  -(—  +  -  +  y  )(x  +  1)  +  ^ 


I  -h  3cos(27rx), 
^  —  3cos(27rx), 


-1  <  X  < 

—  I  <  X  <  0; 
0  <  X  <  |; 

^  <  X  <  1. 


(3.2) 


(3.3) 


This  example  comes  from  [7].  We  use  third,  fifth  and  seventh  order  WENO  schemes  with  101  non- 
uniformly  spaced  points  and  show  the  results  at  t=2  and  8  in  Figure  3.1.  We  can  clearly  observe  better 
resolution  with  increased  order  of  accuracy. 


Example  3.3.  One-dimensional  Burgers  equation: 

{<t>t  +  =0’  -!<»<!, 

I  0(x,  0)  =  —  cos(7rx) 


(3.4) 


with  periodic  boundary  conditions. 

The  local  Lax-Friedrichs  fiux  (2.4)  is  used.  At  t  =  O.S/tt^,  the  solution  is  still  smooth.  We  list  the  errors 
and  the  numerical  orders  of  accuracy  in  Table  3.2,  using  linear  and  WENO  schemes.  We  see  that  the  correct 
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Fig.  3.1.  One- dimensional  linear  equation.  Non-uniform  mesh  with  101  points.  Solid  line  is  the  exact  solution.  Left:  a 
comparison  of  third  order  ‘o  ’  and  fifth  order  ‘A  right:  a  comparison  of  fifth  order  A  ’  and  seventh  order  ‘n  Top:  results  at 
t  — 2;  bottom:  results  at  t  —  %. 


orders  of  accuracy  are  obtained.  At  t  =  3.5/7r^,  the  solution  has  developed  a  discontinuous  derivative.  In 
Figure  3.2,  we  show  the  sharp  corner-like  numerical  solution  with  41  points  using  third,  fifth  and  seventh 
order  WENO  schemes.  From  now  on,  the  solid  line  is  the  exact  solution,  and  the  circles  are  numerical 
solutions. 


Example  3.4.  One-dimensional  equation  with  a  non-convex  fiux: 

Ut  -  cos(<?!>a,  +  1)  =  0,  -l<a:<l, 

I  (p{x,  0)  =  —  cos(7rx) 

with  periodic  boundary  conditions. 

The  local  Lax-Friedrichs  fiux  (2.4)  is  used.  At  t  =  I.S/tt^,  the  solution  has  developed  a  corner-like 
discontinuity  in  the  derivative.  The  numerical  result  with  41  points  is  shown  in  Figure  3.3. 
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Table  3.2 

Accuracy  for  ID  Burgers  equation,  Linear  and  WENO  Sehemes,  t  —  0.5/7r^. 


Linear 

3rd  order 

5th  order 

7th  order 

N 

LOO  qytot 

order 

LOO  QYIOI 

order 

LOO  QYIOI 

order 

10 

8.12E-03 

— 

3.89E-03 

— 

1.75E-03 

— 

20 

1.46E-03 

2.47 

4.29E-04 

3.18 

1.14E-04 

3.95 

40 

2.52E-04 

2.54 

2.33E-05 

4.20 

4.78E-06 

4.57 

80 

3.83E-05 

2.72 

8.94E-07 

4.70 

8.54E-08 

5.81 

160 

5.23E-06 

2.87 

2.94E-08 

4.93 

8.57E-10 

6.64 

320 

6.82E-07 

2.94 

9.39E-10 

4.97 

7.33E-12 

6.87 

WENO 

3rd  order 

5th  order 

7th  order 

N 

LOO  qytoj: 

order 

LOO  qytoj: 

order 

LOO  qyyoy 

order 

10 

5.29E-02 

— 

5.21E-03 

— 

2.92E-03 

— 

20 

1.93E-02 

1.45 

3.78E-04 

3.78 

3.02E-04 

3.27 

40 

4.76E-03 

2.02 

3.05E-05 

3.63 

5.58E-06 

5.76 

80 

5.97E-04 

3.00 

1.26E-06 

4.60 

5.68E-08 

6.62 

160 

3.03E-05 

4.30 

4.19E-08 

4.91 

5.66E-10 

6.65 

320 

l.llE-06 

4.77 

9.18E-10 

5.51 

6.28E-12 

6.50 

Example  3.5.  The  one-dimensional  Riemann  problem  with  a  non-convex  flux: 

Ut  +  i('/>^-l)(0x-4)=O,  -l<a;<l, 

l0(a;,O)  =  -2\x\. 

We  remark  that  this  is  a  demanding  test  case.  Many  schemes  have  poor  resolutions  or  could  even 
converge  to  a  non-viscosity  solution  for  this  case.  The  global  Lax-Friedrichs  flux  (2.4)  is  used.  Local  Lax- 
Friedrichs  flux  is  less  satisfactory  for  this  case.  Numerical  results  at  t  =  1  with  81  non-uniform  grid  points 
are  shown  in  Figure  3.4.  Third,  flfth  and  seventh  order  WENO  schemes  are  used. 


Example  3.6.  Two  dimensional  linear  equation: 


I  +  0a;  +  02/  —  0,  -2  <  X  <  2,  — 2  <  2/  <  2, 

\4>ix,y,0)  =  +  y)). 

with  periodic  boundary  conditions. 

We  first  use  uniform  triangular  meshes,  shown  in  Figure  3.5  for  the  coarsest  case  h  =  ^  where  h  is  the 
length  of  the  right  angled  side,  to  test  the  accuracy  for  the  third  and  fourth  order  linear  schemes  and  WENO 
schemes.  We  then  use  non-uniform  meshes,  shown  in  Figure  3.6  for  the  coarsest  case  N  =  105  where  N  is 
the  total  number  of  nodes  in  the  mesh.  The  refinement  of  the  non-uniform  meshes  is  done  in  a  uniform  way, 
namely  by  cutting  each  triangle  into  four  smaller  similar  ones.  The  accuracy  results  are  shown  only  for  the 
non-uniform  meshes,  to  save  space,  in  Table  3.3.  The  expected  orders  of  accuracy  are  obtained. 
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5th  order,  N=41 


Fig.  3.2.  One- dimensional  Burgers  equation^  loeal  Lax-Friedriehs  flux,  t  =  3.5/7r^.  Solid  line:  the  exaet  solution;  eireles: 
numerieal  solutions.  Non-uniform  mesh  with  points.  Top  left:  third  order  WENO  seheme;  top  right:  fifth  order  WENO 
scheme;  bottom:  seventh  order  WENO  scheme. 


Table  3.3 

Accuracy  for  2D  Linear  equation^  Non-uniform  Meshes ^  Linear  and  WENO  Schemes,  t  — 2. 


3rd  order 

4th  order 

Linear 

WENO 

Linear 

WENO 

N 

j^oo  qytoj: 

order 

LOO  qytoj: 

order 

LOO  qytoj: 

order 

LOO  qyyOY 

order 

105 

3.27E-01 

— 

7.00E-01 

— 

I.IOE-Ol 

— 

6.94E-01 

— 

385 

5.59E-02 

2.55 

2.42E-01 

1.53 

6.99E-03 

3.98 

2.35E-01 

1.56 

1473 

8.04E-03 

2.80 

6.44E-02 

1.91 

4.20E-04 

4.06 

4.72E-02 

2.32 

5761 

1.07E-03 

2.91 

1.16E-02 

2.48 

2.54E-05 

4.05 

3.20E-03 

3.88 

22785 

1.39E-04 

2.94 

5.03E-04 

4.52 

1.49E-06 

4.09 

4.43E-05 

6.18 
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X 

Fig.  3.3.  One  dimension,  non-eonvex,  loeal  Lax-Friedriehs  flux,  H{u)  —  —  cos(tt  +  l),t  =  1.5/7r^.  Solid  line:  the  exaet 
solution;  eireles:  numerieal  solutions.  Non-uniform  mesh  with  points.  Top  left:  third  order  WENO  seheme;  top  right:  fifth 
order  WENO  scheme;  bottom:  seventh  order  WENO  scheme. 


Example  3.7.  Two  dimensional  Burgers  equation: 

{4>‘i  +  4’l  + 1)^ 

^  4>t  +  f ^  =  0, 

</>(»,  y,0)  =  -cos(^!^^t2i). 


-2  <  a;  <  2, -2  <  2/ <  2, 


with  periodic  boundary  conditions. 

At  t  =  O.S/tt^,  the  solution  is  still  smooth.  We  use  both  the  uniform  meshes  (Figure  3.5)  and  the  non- 
uniform  meshes  (Figure  3.6)  to  test  the  accuracy,  but  errors  and  orders  of  accuracy  for  the  third  and  fourth 
order  linear  and  WENO  schemes  only  for  the  non-uniform  meshes  are  listed  in  Table  3.4,  to  save  space. 
We  can  see  that  the  correct  orders  of  accuracy  are  obtained.  At  t  =  the  solution  has  developed 

discontinuous  derivatives.  We  use  the  third  and  fourth  order  linear  and  WENO  schemes  to  compute,  with 
uniform  mesh  ot  h  =  The  results  are  shown  in  Figures  3.7  and  3.8.  We  can  see  that  for  this  example 
which  is  not  very  demanding,  the  linear  schemes  which  do  not  use  WENO  nonlinear  weights  can  also  obtain 
non-oscillatory  results.  It  seems  that  the  nonlinear  WENO  strategy  is  needed  only  for  the  more  demanding 
cases  such  as  those  for  some  non-convex  fluxes. 
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Fig.  3.4.  One- dimensional  Riemann  problem,  global  Lax- Friedrichs  flux,  H{u)  —  \{u‘^  —  —  A),  t  —  1.  Solid  line: 

the  exact  solution;  circles:  numerical  solutions.  Non-uniform  mesh  with  81  points.  Top  left:  third  order  WENO  scheme;  top 
right:  fifth  order  WENO  scheme;  bottom:  seventh  order  WENO  scheme. 


Table  3.4 

Accuracy  for  2D  Burgers  equation.  Non-uniform  Meshes,  Linear  and  WENO  Schemes,  t  —  0.5/7r^. 


3rd  order 

4th  order 

Linear 

WENO 

Linear 

WENO 

N 

order 

order 

order 

LOO 

order 

105 

6.10E-02 

— 

1.68E-01 

— 

2.91E-02 

1.71E-01 

385 

1.68E-02 

1.86 

5.14E-02 

1.71 

6.10E-03 

2.26 

5.61E-02 

1.61 

1473 

3.54E-03 

2.25 

1.51E-02 

1.77 

7.45E-04 

3.03 

1.63E-02 

1.78 

5761 

5.63E-04 

2.65 

2.76E-03 

2.45 

5.66E-05 

3.72 

1.88E-03 

3.12 

22785 

7.70E-05 

2.87 

1.63E-04 

4.08 

2.65E-06 

4.42 

8.64E-05 

4.45 
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Fig.  3.5.  Uniform  mesh  with  h  — 


Example  3.8.  Two-dimensional  equation  with  a  non-convex  flux: 

f  “  cos(0a;  (j)y  1)  =  0,  —2  <  X  <  2,  —2  ^  y  K  2^  (3  9) 

\<p{x,y,0)  =  -cos(^!^^). 

with  periodic  boundary  conditions. 

For  this  example  we  use  a  non-uniform  mesh  (Figure  3.6)  with  N  =  1473  nodes.  At  t  =  the 

solution  develops  a  discontinuous  derivative,  and  we  show  the  results  using  the  third  and  fourth  order  linear 
and  WENO  schemes  in  Figures  3.9  and  3.10.  Similar  to  the  Burgers  equation  case,  linear  schemes  which  do 
not  use  WENO  nonlinear  weights  can  produce  non-oscillatory  results  as  well. 

Example  3.9.  The  two-dimensional  methods  are  applied  to  the  one-dimensional  non-convex  Riemann 
problem  in  Example  3.5.  This  example  is  more  demanding  than  the  previous  two  examples,  and  we  will  see 
that  the  linear  scheme  will  not  work.  We  solve  the  problem  in  Example  3.5  in  the  domain  [—  l,l]x[-0.2,0.2] 
with  the  triangulation  shown  in  Figure  3.11.  The  periodic  boundary  condition  is  applied  in  the  ^/-direction. 
We  plot  the  solutions  along  the  central  cut  line  y  =  0. 

First  we  use  the  third-order  linear  scheme  to  compute  this  problem,  and  reflne  the  mesh  to  test  conver¬ 
gence.  The  results  are  plotted  in  Figure  3.12,  with  h  =  ^5  ^5  We  can  see  that  the  solutions  of  the  linear 
scheme  with  one  mesh  reflnement  does  not  converge  to  the  correct  viscosity  solution. 

Next  we  use  the  third-order  and  fourth-order  WENO  schemes.  From  Figure  3.13,  we  can  see  that 
the  solutions  of  WENO  schemes  converge  to  the  correct  viscosity  solution  when  the  mesh  is  reflned.  And 
obviously,  the  fourth-order  scheme  has  a  better  resolution  than  the  third-order  scheme. 

At  last  we  use  the  first-order  monotone  scheme  with  the  monotone  Lax-Friedrichs  Hamiltonian  (2.2)  to 
compute  the  problem.  The  results  are  shown  in  Figure  3.14.  It  converges  to  the  correct  viscosity  solution. 
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Fig.  3.6.  Non-uniform  mesh  for  the  coarsest  case  with  N  —  105  nodes. 

3rd-order  Linear  Scheme,  h=1  /1 0  3rd-order  WENO  Scheme,  h=1  /1 0 


Fig.  3.7.  Two-dimensional  Burgers  equatiouj  t  =  1.5/7r^.  Uniform  mesh  with  ^  Left:  third  order  linear  scheme; 

right:  third  order  WENO  scheme. 


as  expected.  But  comparing  with  the  results  of  the  third  and  fourth  order  WENO  schemes,  the  first  order 
monotone  scheme  needs  a  much  more  refined  mesh  to  achieve  the  same  resolution. 


Example  3.10.  Two  dimensional  Riemann  problem: 


4>t  +  sin(0a;  +  —  O5 

0(a;,2/,O)  =7T(\y\ -  |a;|). 


-1  <  a;  <  1, -1  <  y  <  1, 


(3.10) 


For  this  example,  we  use  the  uniform  triangular  mesh  with  40  x  40  x  2  elements.  Third  and  fourth  order 
WENO  schemes  are  used.  We  show  the  numerical  results  at  t  =  1  in  Figure  3.15. 
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4th-order  Linear  Scheme,  h=1  /1 0 


4th-order  WENO  Scheme,  h=1/10 


Fig.  3.8.  Two-dimensional  Burgers  equation,  t  =  l.d/n^. 
right:  fourth  order  WENO  scheme. 


Uniform  mesh  with  h  —  Left:  fourth  order  linear  scheme; 


3rd-order  Linear  Scheme,  1473  nodes 


3rd-order  WENO  Scheme,  1473  nodes 


Fig.  3.9.  Two  dimensions,  H{u,v)  —  —cos(u  -\-v-\-l),t  =  1.5/7r^.  Nonuniform  mesh  with  N  =  1473  nodes.  Left:  third 
order  linear  scheme;  right:  third  order  WENO  scheme. 


Example  3.11.  A  problem  from  optimal  control: 


+  (sin^)0a;  +  (sinx  +  ^  sin^  y  —  (1  —  cosx)  =  0,  -tt  <  x  <  tt,  -tt  <  2/  <  tt, 

(t>ix,y,0)  =0. 

(3.11) 


with  periodic  boundary  conditions,  see  [13]. 

We  use  the  uniform  triangle  mesh  with  60  x  60  x  2  elements.  This  means  each  of  the  60  x  60  rectangles 
is  cut  in  two  triangles  diagonally.  Third  and  fourth  order  WENO  schemes  are  used.  The  solution  at  t  =  1 
is  shown  in  Figure  3.16,  and  the  optimal  control  u)  =  sign{(py)  is  shown  in  Figure  3.17. 
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4th-order  Linear  Scheme,  1473  nodes 


Fig.  3.10.  Two  dimensions,  H{u,v)  —  —  cos(tt  +  +  1),  t 
order  linear  scheme;  right:  fourth  order  WENO  scheme. 


4th-order  WENO  Scheme,  1473  nodes 


1.5/7r^.  Nonuniform  mesh  with  N  =  1473  nodes.  Left:  fourth 


Fig.  3.11.  Mesh  for  Example  3.9  with  /i  =  ^. 


Example  3.12.  2D  eikonal  equation  with  a  non-convex  Hamiltonian,  which  arises  in  geometric  optics  [9]: 

=  0<x<l,0<y<l, 

0(x,2/,O)  =  0.25(cos(27rx)  —  l)(cos(27r2/)  —  1)  —  1. 


(3.12) 


We  use  the  third-order  and  fourth  order  WENO  schemes.  We  use  both  the  uniform  meshes  and  the  non- 
uniform  mesh  shown  in  Figure  3.18,  and  present  the  results  in  Figure  3.19  for  the  non-uniform  mesh  case  at 
t  =  0.6. 


Example  3.13.  The  level  set  equation  in  a  domain  with  a  hole: 

<f>t  +  sign{(t>o){yJ 4^1  +  -  1)  =  0,  |  <  ^x'^  +  y‘^  <  1, 

(l>{x,y,0)  =  Mx,y)- 


(3.13) 


This  problem  comes  from  [17].  The  solution  0  to  (3.13)  has  the  same  zero  level  set  as  00?  and  the 
steady  state  solution  is  the  distance  function  to  that  zero  level  curve.  In  this  example,  the  exact  steady  state 
solution  is  the  distance  function  to  the  inner  boundary  of  the  domain.  We  compute  the  time-dependent 
problem  to  reach  a  steady  state  solution,  using  the  exact  solution  for  the  boundary  condition.  The  mesh  is 
shown  in  Figure  3.20,  and  the  result  using  the  fourth-order  WENO  scheme  is  shown  in  Figure  3.21. 
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Fig.  3.12.  Convergence  study  for  Example  3.9,  linear  scheme,  central  line  cut  at  y  —  0.  Solid  line  is  the  exact  solution; 
diamonds  are  for  the  coarsest  mesh  with  h  —  circles  are  for  the  next  refined  mesh  with  h  —  and  reverse  triangles  are 
for  the  most  refined  mesh  with  h—^. 


Example  3.14.  The  problem  of  a  propagating  surface: 

f<^t-(l-eif)y/l  +  02+02  =0,  0<»<1,0<J/<1, 

I  0(^5  2/5  0)  =  1  “  |(cos27rx  —  l)(cos27r2/  —  1), 


(3.14) 


where  K  is  the  mean  curvature  defined  by 

_  ^xxQ-  +  4^y)  —  ‘^^xy^x^y  +  ^yy{^  + 

(1  +  02+02)1  ’ 

and  £  is  a  small  constant.  A  periodic  boundary  condition  is  used. 

This  problem  came  from  [12].  We  use  the  fourth  order  WENO  scheme,  and  the  second  derivative  terms 
are  approximated  by  those  of  the  interpolating  polynomial  using  the  stable  big  stencil  of  our  fourth  order 
scheme  discussed  in  section  2.  For  £  =  0  (pure  convection),  whose  solution  will  develop  a  discontinuous 
derivative  in  the  center  of  the  domain,  we  use  the  non-uniform  mesh  shown  in  Figure  3.18,  and  for  £  =  0.1, 
we  use  the  uniform  mesh  with  50  x  50  x  2  elements  since  the  solution  is  smooth.  The  results  are  shown  in 
Figure  3.22.  Notice  that  the  surfaces  at  t  =  0  and  t  =  0.1  (for  ^  =  0.1)  are  shifted  downward  in  order  to 
show  the  detail  of  the  solution  at  later  time. 


Example  3.15.  A  problem  from  computer  vision: 

{4>t+  I{x,y)^l  +  -  1  =  0,  -1  <  a;  <  1,-1  <  2/ <  1, 

1  </>(«,  2/,0)  =0 


(3.16) 
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Fig.  3.13.  Convergence  study  for  Example  3.9,  WENO  schemes,  central  line  cut  aty  —  t).  Solid  line  is  the  exact  solution; 
diamonds  are  for  the  coarsest  mesh  with  h  =  circles  are  for  the  next  refined  mesh  with  h  =  and  reverse  triangles  are 
for  the  most  refined  mesh  with  h  —  Top:  third  order  WENO  schemes;  bottom:  fourth  order  WENO  schemes. 


where  /  is  defined  by 


1 

yi  +  (1  -  |a;|)2  +  (1  -  12/1)2 


(3.17) 


with  (p  =  0  as  the  boundary  condition. 

This  problem  comes  from  [14].  The  steady  state  solution  of  this  problem  is  the  shape  lighted  by  a  source 
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X 


Fig.  3.14.  Convergence  study  for  Example  3.9,  first  order  monotone  scheme,  central  line  cut  at  y  —  t).  Solid  line  is  the 
exact  solution;  dotted  line  is  for  the  coarsest  mesh  with  h  —  ^  dash-dotted  line  is  for  the  next  refined  mesh  with  h  —  , 

and  the  dashed  line  is  for  the  most  refined  mesh  with  h  —  ^ . 

3rd  order  WENO,  h=1/20  4th  order  WENO,  h=1/20 


Fig.  3.15.  Two-dimensional  Riemann  problem,  H{u,v)  —  sin(w  -\-v),  t  —  1.  Uniform  triangle  mesh  with  h  —  Left: 
third  order  WENO  scheme;  right:  fourth  order  WENO  scheme. 

located  at  infinity  with  vertical  direction.  The  solution  is  not  unique  if  there  are  points  at  which  /(x,  y)  =  1. 
Conditions  must  be  prescribed  at  one  of  those  points  where  I{x,y)  =  1.  The  exact  steady  solution  is 

(/>(a;,j/,oo)  =  (1- |a;|)(l  -  I2/I)  (3.18) 

We  use  the  uniform  mesh  with  50  x  50  x  2  elements  and  the  third-order  and  fourth-order  WENO  schemes 
to  compute  the  steady  state  solution.  The  boundary  conditions  are  imposed  exactly  from  the  above  exact 
steady  solution.  The  results  are  shown  in  Figure  3.23. 
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3rd-order  WENO  scheme,  h=27r/60 


4th-order  WENO  scheme,  h=27r/60 


Fig.  3.16.  Control  problem,  t  —  1.  Uniform  triangle  mesh  with  h  —  Left:  third  order  WENO  scheme;  right:  fourth 
order  WENO  scheme. 


3rd-order  WENO  scheme,  h=27i/60 


4th-order  WENO  scheme,  h=27i/60 


Fig.  3.17.  Control  problem,  t  —  l,uj  —  sign{<py).  Uniform  triangle  mesh  with  h  =  Left:  third  order  WENO  scheme; 
right:  fourth  order  WENO  scheme. 


4.  Concluding  remarks.  We  have  constructed  third  order  and  fourth  order  WENO  schemes  for 
Hamilton- Jacobi  equations  based  on  2D  triangular  meshes.  Extensive  numerical  experiments  have  been 
performed  to  find  the  most  robust  strategies  in  the  choice  and  grouping  of  stencils  for  the  WENO  schemes. 
These  methods  have  high-order  accuracy  for  smooth  problems  and  high-resolution  for  singularity  of  deriva¬ 
tives.  No  parameters  have  to  be  tuned  in  the  algorithms.  Numerical  Examples  are  shown  to  illustrate  the 
capability  of  the  method.  These  methods  should  be  useful  if  geometry  or  adaptivity  requires  unstructured 
meshes  for  the  solution  of  Hamilton- Jacobi  equations. 
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